Exact dynamics and decoherence of two cold bosons in a ID harmonic trap 
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We study dynamics of two interacting ultra cold Bose atoms in a iiarmonic oscillator potential 
in one spatial dimension. Making use of the exact solution of the eigenvalue problem of a particle 
in the delta-like potential we study time evolution of initially separable state of two particles. The 
corresponding time dependent single particle density matrix is obtained and diagonalized and single 
particle orbitals are found. This allows to study decoherence as well as creation of entanglement 
during the dynamics. The evolution of the orbital corresponding to the largest eigenvalue is then 
compared to the evolution according to the Gross-Pitaevskii equation. We show that if initially 
the center of mass and relative degrees of freedom are entangled then the Gross-Pitaevskii equation 
fails to reproduce the exact dynamics and entanglement is produced dynamically. We stress that 
predictions of our study can be verified experimentally in an optical lattice in the low-tunneling 
limit. 

PACS numbers: 03.75Kk, 03.75.Gg, 67.85.-d 
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I. INTRODUCTION 

Theoretical description of a Bose-Einstein condensate 
of trapped weakly interacting atomic system is tradition- 
ally based on a mean field approximation By as- 
suming that many-body wave function can be written 
in a form of N-fold product state, i.e. that all atoms 
occupy the same single particle orbital, the stationary 
Gross-Pitaevskii (GP) equation for the order parame- 
ter is found. Assuming further that the N-fold prod- 
uct approximation holds also in dynamical situations 
one arrives at the time dependent Gross-Pitaevskii equa- 
tion. Under most of experimental conditions there are 
no strong correlations in the system and the GP equa- 
tion turned-out to be extremely fruitful in predicting and 
describing a variety of features of those systems. Soon it 
occurred that also high energy solutions of the GP equa- 
tion can be useful in studying Bose systems at finite tem- 
peratures. The GP equation has become a work horse of 
the theory of weakly interacting ultra cold bosons. 

On the other hand examples when the mean field de- 
scription does not reproduce the real dynamics have been 
studied. For instance direct comparison of the mean field 
and many body theory of vortex nucleation showed that 
the GP equation fails to describe this phenomenon 0- 
0]. Similarly a mean field description of attractive Bose 
systems encounters difficulties p-Q- Due to permanent 
progress in experimental techniques the physics of ultra 
cold atomic gases started to penetrate areas traditionally 
associated with condensed matter physics where correla- 
tions play a crucial role. Evidently, in such situations 
simple mean field description based on the GP equation 
becomes insufficient. The Mott insulator-superfiuid tran- 
sition P or the Tonks- Girardeau gas d, [l^ are some 
examples. 

It is commonly believed that creation of the Mott insu- 



lator with a small and controlled number of atoms per lat- 
tice site allows for applications of such systems in quan- 
tum information. All quantum information processing 
is based on the quantum correlations which cannot be 
described by any classical theory based on local realism. 
States showing these non-local correlations are known as 
entangled states. Two entangled spin-i particles are the 
essence of the Einstein, Podolsky, and Rosen paradox 

The entangled states vialote Bell inequalities 
It has been realized that also continuous variable entan- 
gled systems can be used in quantum computations [l3l - 
[igj . Several authors [20'-^2^ studied recently creation of 
entanglement during dynamics of two interacting parti- 
cles in a harmonic trap. In particular it has been shown 
that this dynamically entangled state violates a Bell-type 
inequality for a certain choice of observables [l^l • 

In this paper we study dynamically created entangle- 
ment, and its measure - the von Neuman entropy, for 
a realistic system of two identical bosonic atoms in a 
harmonic trap. We consider low energy collisions of the 
atoms. At such energies the range of van der Waals in- 
teractions is smaller than the s-wave scattering length. 
Therefore, the interaction potential can be approximated 
by a contact pseudo potential. This approximation oc- 
curred to be in excellent agreement with experimental re- 
sults (23j where binding energy of molecular system has 
been measured. The molecules were created from atoms 
in an optical lattice in the limit of small tunneling. This 
experimental arrangement is perfectly suited for a study 
of exact dynamics of two trapped atoms. We consider 
a realistic case of two atoms per lattice site deep in the 
Mott insulator phase. By applying a Bragg pulses [l^] 
one creates a state in which each atom is in a superposi- 
tion of two counter propagating wave packets. Initially, 
such a state is a two-fold product state of two identical 
wave functions, i.e. is separable. Each wave function 



2 



has two components moving initially with opposite mo- 
menta. The center of mass dynamics is generated by a 
different Hamiltonian than dynamics of the relative co- 
ordinate therefore this two particle continuous variable 
system becomes entangled. We use the von Neuman en- 
tropy as a measure of the entanglement and study its 
behavior in time for various interaction strengths. We 
also analyze a coherence of the system which is directly 
related to the largest eigenvalue of the single particle den- 
sity matrix and compare the exact dynamics to the mean 
field description based on the GP equation. 



II. TWO BOSONS IN A HARMONIC TRAP 

We are going to study dynamics of the simplest non- 
trivial system - two atoms confined in a one dimensional 
harmonic potential. In fact generalization of our results 
to two or three spatial dimensions is straightforward. We 
limit our analysis to the ID case as this situation cap- 
tures all features of the dynamics. For simplicity we are 
using harmonic-oscillator units. It means that all ener- 
gies are measured in hcj, all lengths in -y/ h/mut, and all 
momenta in yhmuj. Hamiltonian of the system of two 
interacting bosons in the harmonic trap has the form: 
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where xi and X2 are positions of atoms interacting via 
a short range potential modeled by the delta function. 
This form of the short range interaction is justified in 
the limit of vanishing relative velocity of colliding atoms, 
where atomic de Broglie wavelength is much larger then 
a range of two body potential. In 2D and 3D the corre- 
sponding Hamiltonian is not a self adjoint operator. To 
correct for this fact a regularization is required. In con- 
trast to many dimensions, the regularization of the delta 
function is not necessary in one dimensional case [2^ . In 
ID the parameter g is given hy g = — 2/ao, where ao 
is a scattering length ,2(|]. It is worth to notice that fi- 
nite range interactions between particles modeled by the 
Gaussian function in the context of the dynamics of two 
bosons was studied in [27| . 

To demonstrate entanglement formation we study the 
evolution of two bosons which initially are in a product 
quantum state 



*o(a;i,a;2) = *o(a;i)$o(a;2)- 



(2) 



center of mass and the relative coordinates: 
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(3a) 
(3b) 



In these coordinates Hamiltonian ([T} separates into two 
independent parts - the center of mass part, and the 
relative part: 
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As one can see, the dynamics of the center of mass is 
described by the standard one dimensional harmonic os- 
cillator Hamiltonian ([^a]). Its eigenstates are well known 
and have a standard form 



X„(X) = ^H„(X)e-^V2^ 
V2"n! 



(5a) 



where H„(x) are Hermite polynomials. The energy of 
n-th eigenstate in our units is obviously given by 
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(5b) 



The eigenstates of the Hamiltonian (j4bl) describing rela- 
tive dynamics of two particles are given in ^28i] and for 
one dimensional problem have a form 
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Hm(C)e~^'/^ TO odd (6a) 



^rniO = U(-z^,„, C^) & TO cven (6b) 

where U(a, /?, x) are confluent hypergeometric functions, 
and Affn are normalization coefficients. Since the wave 
function of identical bosons must be symmetric under ex- 
change of the two particles, therefore the physical wave 
function is composed from functions with even m only. 
The energies Em of these even states are given by a se- 
quence of zeros of the function: 



fiE) = 



r(-£;/2 + 3/4) 1 
r(-£;/2 + i/4) " 



(7) 



The quantum number v is equal to Vm = (2i?m — l)/4. 
The initial wave function can be easily decomposed to 
the superposition of the eigenstates of the Hamiltonian: 



Function $0(2;) is a one-particle wave function called the 
order parameter in the mean field context. 

The exact dynamics of the two interacting bosons in 
the harmonic trap can be found because all eigenstates 
of the full two-body Hamiltonian ([l} are known. They 
are found in |28|. The two particle problem has to be 
first brought to a single particle one by introducing the 
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Obviously the evolution of the initial two boson state is 
given by: 



v|/(e,X,t) = V 
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(9) 
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The last step is to return to the original coordinates by 
using relations ([3]): 



X ^„ (^i^j ^-^iSMt (10) 

Standard method of detection of ultra cold trapped 
atomic systems are destructive. The optical lattice po- 
tential is turned off and the system is allowed to expand 
balistically. Only after expansion a size of the system 
exceeds a resolution of a CCD camera. The picture of 
the CCD camera gives therefore direct insight into the 
initial momentum distribution of atoms. The wave func- 
tion Eq. ljlOl) written in the momentum space of the two 
atoms is: 

/■oo POO 



(11) 

In repeated single particle detections preceded by the 
ballistic expansion of the system one-particle momentum 
distribution is monitored; 



?^Exact(fc,^) = pi{k, k,t), 



(12) 



where pi{k, k',t) is the reduced one particle density ma- 
trix in the momentum representation: 



pi{k,k',t) 



dk2ip*{k,k2,t)'ilj{k',k2,t) (13) 



By making its spectral decomposition we can determine 
the number of orbitals and their relative occupations 
needed for accurate description of the two bosons dy- 
namics. Time dependence of the eigenvalues of the den- 
sity matrix is discussed below. Let us mention that the 
largest eigenvalue is a direct measure of the coherence of 
the system. 

We shall compare this exact dynamics with the ap- 
proximate one governed by the Gross-Pitaevskii equa- 
tion. The main idea leading to the mean field approxi- 
mation relies on the assumption that generation of entan- 
glement between bosons during the evolution is negligible 
and therefore the quantum state of the system remains 
separable. In other words all correlations between bosons 
are neglected and the same wave functions of every par- 
ticle is assumed during the evolution: 



*(xi,X2,t) = $(xi,i)$(x2,t). 



(14) 



This assumption leads directly to the Gross-Pitaevskii 
equation which determines the dynamics of the one- 
particle wave function $(x,t): 



idMx,t) = ( + ^x^+gMxA)\^ ] <i>{x,t). 



The probability density in momentum space reads: 

nMk,t) = mk,t)\\ (16) 

where (f>{k,t) is the Fourier transform of the time 
dependent solution of the CP equation, 4'{k,t) = 
J da; e~''^^$(a;, i). We compare the exact one-particle 
momentum distribution with that predicted by the 
Gross-Pitaevskii approximation ([T6l). Moreover, in the 
situation when many eigenvalues of the density matrix 
are of the same order we can also compare the Gross- 
Pitaevskii momentum distribution (|16p with the momen- 
tum distribution of the dominant orbital obtained from 
diagonalization of the exact one-particle density matrix 
in the momentum space. Obviously, in a general case the 
GP solution overestimates the coherence of the system. 

The GP equation is solved numerically on a spatial grid 
of Np = 2^° points separated by Sx = 5 ■ 10~^. The time 
step is equal to St = 7r/4-10~^. We use the split-operator 
method which is very stable for the chosen temporal and 
spatial steps. 



III. RESULTS 

To make the detailed comparison we concentrate on a 
one particular class of the initial states. We assume that 
at the beginning each particle is in the state described 
by the Schrodinger cat like wave function 



(17) 



(15) 



Parameter L measures the separation between two wave 
packets moving in the opposite direction in the relative 
coordinates space. Such a choice is motivated by the 
preparation procedure described above, i.e. we assume 
that Bragg pulses bringing the atoms into the superpo- 
sition of wave packets moving in opposite directions are 
applied. When L — then the initial state is very close 
to the ground state of the system so we expect that the 
exact dynamics is almost indistinguishable from the dy- 
namics in the mean filed approximation. For large L the 
initial state is still separable but it is highly delocalized. 
Relative and center of mass degrees of freedom are en- 
tangled in the initial state. They evolve in a different 
way, therefore we expect that the exact dynamics could 
be dramatically different than the dynamics predicted by 
a simple mean field approach. 



A. Dependence on delocalization of one-particle 
state 

First let us discuss situation for generic interaction 
strength g = —0.2 (a = 10) when L — 1, i.e. when the 
extension of the initial state is equal to the trap length 
unit. We observe that the single particle density matrix 
obtained from the exact dynamics develops more then 
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FIG. 1: (a) Eigenvalues of the one-particle density matrix 
(|13|l . Unit of time is equal to the period of the trap. In this 
situation (parameters: g = —0.2, L = 1) the initial state 
is not far from the ground state of the system. One eigen- 
value still dominates, therefore system should be quite well 
described by the mean field approximation, (b) Two plots 
present the one-particle momentum distributions predicted 
by the exact (thick solid line) and the Gross-Pitaevskii solu- 
tions (dotted line) in two interesting moments. Third (thin 
solid) line comes from the exact solution and presents the mo- 
mentum distribution of the first orbital. As was expected all 
three predictions are almost the same for considered set of 
parameters. Movie presenting time evolution of momentum 
distributions is available on-line [291 



one nonzero eigenvalue, i.e. many one particle orbitals 
are involved. Fig. [T^ shows time dependence of the eigen- 
values of the one-particle density matrix (IT51) . Because 
one of the eigenvalues is incessantly much larger than 
the others the system coherence is large and the Gross- 
Pitaevskii description is quite accurate in this case. Time 
dependence of the momentum distributions deduced from 
the Gross-Pitaevskii equation and the exact solution are 
shown in Fig. [T}3 and they are in agreement with our 
predictions (whole time dependence of momentum dis- 
tributions is available on-line [1^). 

Situation changes dramatically when we increase the 
delocalization parameter. When L is large enough then 
a few orbitals can play the crucial role in the dynam- 
ics and the mean field approximation is no longer valid. 
Fig. [5^ shows the time dependence of the eigenvalues of 
the density matrix for L = 3. As we see, the main or- 
bital (its eigenvalue is represented by a thick solid line) 
initially dominates. But after a few periods of the trap 
oscillations the other orbital becomes much more impor- 
tant than the first one. The dynamics is obviously much 
more complicated than it is predicted by the mean field 



approach. It is clear when we compare the momentum 
density distribution predicted by the exact and the mean 
field solutions (Fig. ^Bp and movie available on-line poj). 

We see that evidently Gross-Pitaevskii equation prop- 
erly describes the dynamics of the first orbital rather then 
the whole system. Fig. [^b- It is the reason why the Gross- 
Pitaevskii equation gives good predictions when only one 
eigenvalue of the one particle matrix dominates during 
the entire evolution. 

(a) , , 
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FIG. 2: (a) Eigenvalues of the one-particle density matrix (|13p 
for g = —0.2, L = 3. Unit of time is equal to the period of the 
trap. In this situation the initial state is a product of highly 
delocalized one-particle wave functions. There is no one dom- 
inant eigenvalue during the evolution and therefore the Gross- 
Pitaevskii equation will not predict dynamics correctly, (b) 
Time dependence of the one-particle momentum distribution 
predicted by the exact (thick solid line) and Gross-Pitaevskii 
solutions (doted line). As long as the first eigenvalue domi- 
nates during the time evolution the predictions are almost the 
same. After five periods the second eigenvalue is the largest 
one and therefore the predictions are highly different. So- 
lutions of the exact and GP dynamics become similar after 
eleven trap periods when the first eigenvalue starts to domi- 
nate again. Notice that third (thin solid) line presenting mo- 
mentum distribution of first one-particle orbital of an exact 
solution recovers predictions of the Gross-Pitaevskii equation. 
The movie is available on-line |29| . 

It is also interesting to study similarities and differ- 
ences between prediction of the mean field approach and 
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FIG. 3: Time dependence of the one-particle momentum 
distribution for the antisymmetric initial state ^oix) — 
-(a,-L)V2 _g-(^+L)2/2j ^j^jj L = 3. Thick sohd line 

represents the density predicted by an exact solution, while 
doted one the density coming from the mean field approach. 
Properties of the Gross-Pitaevskii equation provide that if the 
function $o{x) is antisymmetric under x — >■ —x symmetry 
then it will stay antisymmetric during the whole evolution. It 
is not true for the one-particle density predicted by an exact 
solution. Thin solid line presents momentum distribution of 
the first one-particle orbital of the exact solution. As we see 
its spatial reversal symmetry is preserved during the evolu- 
tion. It shows once more that Gross-Pitaevskii equation de- 
scribes properly the dynamics of the first orbital only. Whole 
movie is available on-line i29i]. 



exact solution in the situation when the initial state of 
each particle is antisymmetric in position space, i.e. is 
described by the wave function of the form 



-{x+Lf/2 



(18) 



Nevertheless the wave function of the system is still sym- 
metric under particle exchange. Since corresponding GP 
Hamiltonian is invariant under reflection x — >■ — a;, there- 
fore the symmetry of the initial state will be preserved. 
As we observe, it is not true for the exact two body 
dynamics. The evolution preserves only the symmetry 
of each orbital separately, but not the symmetry of the 
whole system. It is clearly demonstrated in Fig. [3] where 
we compare the momentum distribution predicted by the 
mean field approach with the single particle density ob- 
tained from the exact dynamics. Time dependence of the 
eigenvalues of one-particle density matrix is identical as 
the one for the corresponding symmetric case (Fig. [5^). 
Whole movie is also available on-line 129i . 
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FIG. 4: Eigenvalues of the one-particle density matrix as func- 
tions of time for parameters: g = —0.4, L = 2. The inter- 
action between bosons is strong and the initial state of one 
particle is highly delocalized. In such a situation many or- 
bitals play a crucial role during the evolution of the system. 
Therefore the exact dynamics is much more complicated that 
the dynamics predicted by the mean field approximation. 



the situation when the interaction is switched off, the 
two particles initially in the state which is not entangled 
(product state) will stay in such a state during the whole 
evolution even for a highly delocalized state. In this case 
the mean field approximation naturally leads to the same 
solution as the exact solution. It is the interparticle in- 
teraction which can produce entangled two body states 
during the evolution. 

Time dependence of the eigenvalues for a moderate in- 
teraction strength {g = —0.2) is presented in Fig. [T]and 
[21 In those situations only two eigenvalues (i.e. two or- 
bitals) are important for many trap periods. For stronger 
interactions this picture changes significantly. Time evo- 
lution of eigenvalues for strong interaction {g = —0.4) 
and L = 2 are shown in Fig. 21 After a few trap peri- 
ods many different orbitals become important. Moreover 
the orbital which dominates at the beginning becomes 
unimportant after a very short time. Therefore we do 
not expect that the Gross-Pitaevskii approximation may 
give correct predictions in this case. 

On the other hand when the interaction is very weak 
we can expect that the production of entanglement will 
be very slow even for highly delocalized states and there- 
fore the mean field approximation may be correct for a 
long evolution time. Time dependence of the eigenvalues 
of the one-particle density matrix when the interaction 
is weak but the initial state is highly delocalized is pre- 
sented in fig. O 



B. Dependence on interaction strength 

Now we want to show that correctness of the mean 
field approximation significantly depends on the interac- 
tion strength parameter g. It is quite obvious that in 



C. Revivals of product states 

Looking at fig. [51 and fig. [SI one can observe that as 
initially only one eigenvalue dominates in the Schmidt 
decomposition of the single particle density matrix, the 
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FIG. 5: Eigenvalues of the one-particle density matrix as func- 
tions of time for parameters: g = —0.04, L — 2. In this sit- 
uation the interaction between bosons is very week but the 
initial state is far from the ground state of the system. During 
the first eighteen trap periods only one eigenvalue dominates, 
therefore the dynamics of the system can be quite correctly 
described by the mean field approximation for a long time. 
Notice that time scale is two times larger than in the previ- 
ous situation. 



other eigenvalues become more important at later times. 
However, the time dependence of the dominant eigen- 
value exhibits some oscillations and a partial revival of 
the 'initial' eigenvalue can be observed. It is interesting 
to find a physical explanation of this behavior. To this 
end in fig. [6] we show the spectrum of the two-particle 
state for the two studied parameter sets g = —0.2 and 
L = 3 (as in fig. [5]) and g = —0.4 and L = 2 (as in fig. O. 
In this figure we plot the probability of the given eigenen- 
ergy, |a„mp, resulting from the decomposition dH). 

First let us notice that the eigenenergies do appear 'in 
pairs'. The effect of 'pairing' of eigenenergies can be eas- 
ily explained. For given n and m the eigenenergy has 
two components. is the energy of the center of mass 
while Em corresponds to the relative coordinate. As was 
mentioned before, energies of the relative excitations are 
very close to the energies of harmonic trap. This is be- 
cause the potential in the relative coordinate space is the 
harmonic potential of the trap modified at ^ = by the 
presence of delta function. This delta function shifts the 
harmonic energy by very small amount. Therefore state 
labeled by (n, m) is almost degenerate with the state la- 
beled by (m, n) and therefore spectrum is paired. 

If by A we denote the difference between two energies 
of dominant pair (maximum in the spectral decomposi- 
tion) it is clear that after time Tr = 27r/A these two the 
most important amplitudes of the two-atom wavefunc- 
tion match in phase and the partial revival of the product 
state can be observed. This is signified by a reappearance 
of the initial eigenvalue of the single particle matrix. The 
revival time calculated this way for parameters of Fig. [5] 
is equal to Tr = 10.96 while for Fig. His Tr = 5.56 which 
agrees perfectly with predictions of the exact dynamics. 

Obviously revival time tj^ depends on the initial state 



FIG. 6: Energy spectrum of two initial states described by 
the parameters: g = —0.4 and L = 2 (crosses); g = —0.2 and 
L = 3 (squares). Note that both spectra are well picked and 
energies are 'paired'. 

as well as eigenmodes of the Hamiltionian ([1]) . 

In addition to this large time oscillations of the eigen- 
values of the density matrix some small fast oscillations 
can be also observed. This fast modulations appear ev- 
ery half of the trap period when two wavepackets meet 
at the trap center and results from interaction between 
them. 



D. Entanglement of particles 

Mutual interactions between particles obviously leads 
to the quantum correlations between particles. To study 
them we use the correlation measure introduced in j3C| : 

K(pi)-(^E^') (19) 

where are the eigenvalues of the one-particle density 
matrix p\ . This measure has very simple interpretation. 
It gives an effective number of single particle orbitals oc- 
cupied in the given many body state. In particular when 
one-particle density matrix has n equal eigenvalues then 
K = n. 

Other commonly used [20l - [2^ IsTj measure of entan- 
glement in the system is von Neumann entropy defined 
as 

S(pi) = -Tr (pilog pi) = - ^ A, log A, (20) 

i 

This entropy is even more interesting than the number 
of dominant eigenvalues K since it is directly connected 
with the entropy defined in thermodynamical context. 
Time dependence of this two measures of entanglement 
in the system for two different regimes of interaction 
strength are presented in Fig. [71 Obviously in the begin- 
ning, when the system is in separable state, entanglement 
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FIG. 7; Time dependence of the number of dominant eigen- 
values K defined in (|19|) and of the von Neumann entropy S 
defined in ((20]) for g = -0.04 (thick line) and g = -0.2 (thin 
line). Other parameters are the same as in Fig. [S] 



and von Neumann entropy are equal to 1 and respec- 
tively. We observe that correlation K and entropy S in- 
crease in time and seem to saturate for large time. Even 
though they have different physical interpretation they 
behave very similarly which might seem quite surprising. 
They reach 'stationary regime' faster for stronger inter- 
actions. 

Both quantities exhibit fast oscillations modulated by 
a slowly varying functions. These fast oscillations are re- 
lated to partial revival of dominant eigenvalue discussed 
in previous subsection. Every minimum observed in cor- 
relation function corresponds to the moment when there 
is a dominant eigenvalue in the Schmidt decomposition of 
the one-particle density matrix. Let us remind that this 
revivals are related to phase matching of two dominant 
eigenmodes of the two particle state. 

Long time modulations of correlation functions are re- 
lated to the quantum nature of the system and discreet- 
ness of the energy spectrum. In such cases evolution is 
always quasi-periodic and due to the interference of am- 
plitudes long time scale oscillations do appear. In our 
case the number of modes with no zero amplitudes is 
relatively small and therefore oscillations of correlation 
functions appear on a time scale of few hundred trap pe- 
riods. 



contact potential. We assumed that initially each parti- 
cle is transferred by the Bragg pulses to the state being 
the superposition of two wave packets moving in opposite 
directions. We show that the two particle state, although 
initially being a product state does not preserve the prod- 
uct form during the evolution. The reason is that the 
initial state entangles the center of mass and relative co- 
ordinates of the two particle system. These two degrees 
of freedom evolve according to different Hamiltonians. 
As a result the single particle reduced matrix develops 
many eigenvalues during the evolution what signifies de- 
creasing coherence of the system. This situation cannot 
be correctly described by the GP equation. Our predic- 
tions can be verified in the experiment with deep optical 
lattices when two atoms occupy each site. We show one- 
particle momenta distributions for different initial states 
and compare them to those obtained from the mean-field 
dynamics. The differences between the two signify the 
two atom entanglement. The momentum distribution is 
directly measured by exposure of the system to a reso- 
nant light after ballistic expansion and therefore creation 
of entanglement in the two particle system can be eas- 
ily traced in time and compared to exact solutions. We 
monitored the von Neumann entropy which is common 
measure of entanglement. We show that entanglement 
is dynamically created during evolution, however it is 
not very surprising for interacting system. A comment 
about a system of two fermionic particles would be also 
in place. As two identical fermions do not interact in 
the s-wave channel, as long as other partial waves can be 
neglected, their dynamic is driven by the noninteracting 
Hamiltonian. On the other hand, if the spatial part of 
a wave function of two fermions is symmetric and a spin 
part is responsible for the antisymmetrization of the total 
wave function, then our exact solution evidently applies 
to such a situation. 
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